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Abstract. In this review, we will discuss the results of our recent work [T] to study the general 
optimization of the pure isovector parameters of the popular relativistic mean-field (RMF) and 
Skyrme-Hartree-Fock (SHF) nuclear energy-density functionals (EDFs), using constraints on the 
pure neutron matter (PNM) equation of state (EoS) from recent ab initio calculations. By using 
RMF and SHF parameterizations that give equivalent predictions for ground-state properties 
of doubly magic nuclei and properties of symmetric nuclear matter (SNM) and PNM, we found 
that such optimization leads to broadly consistent symmetry energy J and its slope parameter L 
at saturation density within a tight range of a (J) < 2 MeV and a(L) < 6 MeV. We demonstrate 
that a clear model dependence shows up (a) in the curvature parameter of the symmetry energy 
-Ksym, (b) the symmetry energy at supra-saturation densities, and (c) the radius of neutron 
stars. 



1. Introduction 

Phenomenological nuclear effective interactions offer a compact description of the in-medium 
nucleon-nucleon interaction and are useful tools in the applications of both the nuclear structure 
and the astrophysical phenomena. The effective interaction is typically dependent on few 
parameters representing, for example, coupling constants, which are often fit to well-determined 
experimental nuclear observables such as binding energies, charge radii, single particle energy 
spectra and spectra of collective excitations. One of the main objective of modern nuclear many- 
body theory is to obtain an EDF [2] with clear physical connections to ab initio nucleon-nucleon 
interactions and QCD. 

In the recent years, much effort has been devoted to constrain the energy per neutron of PNM 
(.Epnm) at sub-saturation densities. By studying the universal behavior of resonant Fermi gases 
with infinite scattering length, a significant constraint is achieved for the EoS of dilute neutron 
matter [3]. These calculations have been extended to higher densities using the full power 
of quantum Monte Carlo methods [H [5]. Moreover, by studying the physics of chiral three- 
nucleon forces the EoS of PNM is obtained perturbatively up to nuclear saturation density |6] . 
Finally, the auxiliary field diffusion Monte Carlo (AFDMC) technique, which takes into account 
the realistic nuclear Hamiltonian containing modern two- and three-body interactions of the 



Argonne potential and Urbana family of three-body nucleon forces, is used to calculate the EoS 
of PNM up to and above saturation density El [9] . 

In this work we concentrate on the widely used RMF |10(, 111] and SHF |1 2(, 1 13j models, with the 
latter thought of as a non-relativistic expansion of the former |14|I15|. Both models have less than 
ten free parameters in their simplest forms. Although there are more than 200 parameter sets in 
the literature for the SHF model [16] and dozens of parameterizations exist in the simplest form 
of the RMF model, e.g. [T7], many of these sets are old and have been superseded by parameter 
sets fit to more recent and accurate data. Since the number of experimental observables is usually 
larger than the number of free parameters, the problem of optimizing these EDFs is generally 
overdetermined, and this results in a significant degeneracy among parameter sets. Fortunately, 
one can use the covariance analysis techniques [THJ, Q15] to study correlations between predicted 
observables from a particular EDF in its model space. We use the linear regression method 
to optimize the two pure isovector parameters of RMF and SHF models by using the results 
from the ab initio theoretical calculations of the PNM EoS as our 'experimental' constraints. To 
obtain meaningful theoretical uncertainties for the model parameters, as well as for the predicted 
observables, we employ the covariance analysis technique. Optimizing the two pure isovector 
parameters will ensure that the predictions for the well-determined isoscalar observables such as 
binding energies B(A) and charge radii R c h of doubly magic nuclei will not be affected. We find 
the best fit values and la confidence intervals on the properties of isospin-asymmetric nuclear 
matter, such as the symmetry energy parameters. In addition, poorly constrained observables 
such as the neutron skin thickness of lead and neutron-star radii are predicted from the resulting 
constraints. We discuss the manifestation of the model dependence in our results, by exploring 
the symmetry energy at supra-saturation densities and properties of neutron stars. We should 
note that our aim is not to establish new parameterizations of these EDFs, or to set absolute 
constraints on symmetry energy, but to explore as far as possible the generic constraints that 
can be placed by each model on neutron-rich systems once constrained by the information from 
the PNM EoS. 

2. Nuclear Symmetry Energy 

The nuclear symmetry energy, S(p), is defined as the coefficient of the leading term of isospin 
asymmetry parameter, a = (p n — p p )/p, in the expression of the binding energy per nucleon in 
neutron-rich nuclear matter 

E(p, a) = E (p) + S(p)a 2 + 0(a 4 ) , (1) 

where p is the baryon number density with p n (p p ) being the neutron (proton) number density. 
Around the saturation density po, one can express the symmetry energy as 

S(p) = J + L X +^K symX 2 + 0(x 3 ) , (2) 

where \= (p — po) /3po, J is the value of the symmetry energy at saturation density, L is the 
slope parameter, and K syra is the curvature parameter of the symmetry energy at saturation 
density. The coefficients of the higher-order terms in Eq. (pQ) are generally much smaller than 
S(p), so it is usually a good approximation to write the symmetry energy as the difference 
between the energy per nucleon of PNM and SNM, i.e., S(p) ~ E(p, 1) — Eq{p). However, in 
this work we will not use such an approximation, but rather calculate it from the full analytical 
expression in a given model. 

3. Relativistic Mean-Field and Skyrme-Hartree-Fock Models 

We apply the constraints on the microscopic PNM calculations to the two popular 
phenomenological nuclear many-body models to study the nuclear symmetry energy and related 



quantities of nuclear physics and nuclear astrophysics. The commonly used RMF model contains 
an isodoublet nucleon field (?/>) interacting via the exchange of the scalar-isoscalar cr-meson 
((f)), the vector-isoscalar w-meson (V**), the vector-isovector p- meson (b^), and the photon 
(^4^) [10 ] 1114 [20] . The effective Lagrangian density for the model can be written as 
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d fl A u —d u A^ l are the isoscalar, isovector, 
and electromagnetic field tensors, respectively. The nucleon mass M and meson masses m s , m v , 
and m p may be treated (if wished) as empirical parameters. The effective potential U((j>, V^,b M ) 
consists of non-linear meson interactions that simulates the complicated dynamics encoded in few 
model parameters. In the present work we use the following form of the effective potential 



U^V^b^) = -(g s 



(4) 



This model is described by seven interaction parameters: {g s ,g v ,9p, K, A, £, A v }. 

The standard form of the energy density obtained from the zero-range Skyrme interaction 
using the Hartree-Fock method can be written as [22] 
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Here p q , T q , and J g (g = p,n) are, respectively, the number, kinetic energy, and spin-current 
densities, and p, r, and J are the corresponding total densities. The SHF model is expressed 
in terms of nine Skyrme parameters: {to,ti,t 2 ,t3, xo, xi, x 2 , x^, a} and the spin-orbit coupling 
constant Wq, which is taken as 133.3 MeVfm 5 [23] in the present work. 



4. Linear Regression and Covariance Analysis Method 

The linear regression and covariance analysis method is discussed in details in Ref. [25] and its 
power has been recently illustrated in Refs. [HI [191 El]- In a nutshell, one can describe it as 
follows. First, an optimal parameters set is found for a given phenomenological model On h \p) 
as a function of the F model parameters p = (pi, ■ ■ ■ ,Pf)i by minimizing the quality measure x 2 
through a method of a least-squares fit: 

* w= £l — m. — J • (6) 

where N is the number of experimental observables oi exp ^ that are determined with an accuracy 
of AO n . Then one can compute the symmetric matrix of second derivatives: 



x\p) ~ X^Po) = A X 2 (x) = x T A4 x , 



(7) 



where x; = iE_ Ssli. The matrix jCi contains all the information about the behavior of the 

1 (Po)i 

X 2 function around the minimum. In particular, one can find the meaningful theoretical 
uncertainties by computing the statistical covariance of two observables A and B defined as 
follows: 

1 <«> 

The variance & 2 (A) of a given observable A is then simply given by a 2 (A) = cov(A, A). Finally, 
one can plot the covariance ellipses between two observables A and B by diagonalizing the 2x2 
covariance matrix: 

a_( cov{A,A) cov{A,B) \ 

^ cov(B,A) cov(B,B) ) 1 > 

The eigenvalues of this covariance matrix represent the semi-major and semi-minor axes of the 
covariance ellipse, while the eigenvectors provide the orientation of the ellipse. 

For our set of 'experimental' observables oi exp ^ in the \ 2 input we choose the theoretical 
microscopic calculations of the energy per neutron £?pnm m the density range of 0.04 < p < 0.16 
fm~ 3 [B1E1 ES| . We restrict our experimental input to this band only, where the upper bound 
is the PNM results at saturation density, as the extension of the PNM calculations to higher 
densities using piecewise polytropes in the chiral effective theory [27] was shown to allow a 
huge uncertainty window in the EoS. Moreover, the symmetry energy coefficients should only 
be sensitive to the EoS around the saturation density. 



5. Results 

We first identify that the two parameters in each model — g p and A v in the RMF model [28J , and 
xq and £3 in the SHF model — are solely isovector parameters. The change of these parameters 
affects only the isovector sensitive properties of nuclear matter, such as the symmetry energy 
S(p), while the EoS of SNM and therefore properties of SNM, such as saturation density po, 
binding energy per nucleon at saturation density Eq, incompressibility coefficient at saturation 
density Kq, and isoscalar effective mass M* at saturation all remain unchanged. Therefore we 
optimize these two isovector parameters [F = 2 in Eqs. ([7]) and ( [8]) ] with respect to the available 
range of PNM equations of state to constrain the values of the symmetry energy parameters 
at saturation density by using the linear regression and covariance analysis method discussed 
in the previous section. Note that most properties of the SNM are constrained experimentally 
within less than 10%. However, the EoS of PNM predicted by various parameterizations of 
both models differ significantly. Although some of them fall within the band of the microscopic 
PNM calculations, for most of the parameterizations there is almost no or very little agreement 
with these calculations as can be seen from Fig. 1 (a). Therefore, in general, the RMF and the 
SHF model predictions for the symmetry energy parameters are substantially different. Once 
the two isovector parameters are optimized to the energy per neutron .Epnm predictions at sub- 
saturation densities [See Fig. 1 (b)], we find that the symmetry energy parameters at saturation 
can also be significantly constrained (See Table [p. 

To further assess theoretical uncertainties in the symmetry energy parameters, we select 
the accurately-calibrated NL3* [29] and the recent IU-FSU [30] parametrizations from the RMF 
model. The IU-FSU is the recent parametrization that predicts a very soft symmetry energy, and 
was validated against experimental, observational, and theoretical data, while the accurately- 
calibrated NL3* parametrization gives a much stiffer EoS of both SNM (larger value of Kq and 
smaller value of £ parameter) and a stiff symmetry energy (larger values of symmetry energy J 
and slope L) and therefore offers a suitable contrast to IU-FSU. 

To compare the RMF and SHF models on the same basis, we create two SHF parametrizations 
which give the same properties of nuclear matter at saturation as the two RMF parametrizations 
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Figure 1. Comparing the PNM EoS from 11 RMF (blue lines) and 73 SHF parameterizations 
(red lines) that are produced since 1995, with the AFDMC EoS in the AV8'+UIX 
Hamiltonian the variational APR EoS [26J, the low-density band from the constraints 
of resonant Fermi gases [3], and the high-density band from the chiral effective field theory 
calculations with 3- neutron forces [5], before (a) and after (b) the PNM optimization. 



Table 1. Predicted ranges for symmetry energy parameters within RMF and SHF models 
before (with superscript '0') and after (without superscript '0') their pure isovector parameters 
optimized to PNM, and taking into account all remaining variation from parameterizations 
constructed since 1995. 





J u (MeV) 


J (MeV) 


L u (MeV) 


L (MeV) 


(MeV) 


K T (MeV) 


RMF 
SHF 


30.3 - 38.7 
27.8 - 39.6 


30.2 - 31.4 
30.1 - 33.2 


47.2 - 122.7 
5.8 - 100.1 


36.1 - 59.3 
28.5 - 64.4 


-701.7 - -195.3 
-514.8 - -266.3 


-329.7 - -215.7 
-418.8 - -235.3 



(See Table [2]), through the method of writing the Skyrme parameters as functions of macroscopic 
nuclear quantities |23[ [31] . These new parametrizations are referred to as SkNL3* and SklU-FSU 
forces pQ. 

Several definitions of the nucleon effective mass exist in the literature [32]. In the RMF 
model the Dirac effective mass is defined through the scalar part of the nucleon self-energy in 
the Dirac equation. It has been well documented that there is a strong correlation between the 
Dirac effective nucleon mass at saturation density and the strength of the spin-orbit force 
in nuclei [T4"l l33j 134] . Indeed, one of the most compelling features of RMF models is the 



Table 2. Macroscopic quantities from four reference parameterizations. They are the nuclear 
saturation density po, the binding energy per nucleon Eq and incompressibility Kq of SNM at 
saturation, the symmetry energy J and its slope parameter L at saturation, and the nucleon 
effective mass M* at saturation. For consistency, we present the Lorentz effective mass of the 
RMF model, which is set equal to the isovector and isoscalar effective masses in the SHF model. 





Po (fm" a ) 


E (MeV) 


K (MeV) 


J (MeV) 


L (MeV) 


M* (M) 


NL3* 


0.1500 


-16.32 


258.49 


38.7 


122.7 


0.671 


SkNL3* 


0.1527 


-15.76 


258.49 


38.7 


122.7 


0.671 


IU-FSU 


0.1546 


-16.40 


231.33 


31.3 


47.2 


0.687 


SklU-FSU 


0.1575 


-15.70 


231.33 


31.3 


47.2 


0.687 



reproduction of the spin-orbit splittings in finite nuclei. It is shown that models with effective 
masses outside the range 0.58 < MX/M < 0.64 will not be able to reproduce empirical spin- 
orbit couplings [35], when no tensor couplings are taken into account. On the other hand, the 
non-relativistic effective mass parameterizes the momentum dependence of the single particle 
potential, which is the result of a quadratic parametrization of the single particle spectrum. 
A recent study [IB] puts a bound of 0.69 < M*/M < 1.0 for the non-relativistic effective 
masses. It has been argued |36j that the so-called Lorentz mass M£ should be compared 
with the non-relativistic effective mass extracted from analyses carried out in the framework 
of nonrelativistic optical and shell models. For consistency, we choose the effective mass in the 
SHF parameterizations to be equal to the Lorenz mass in the RMF parameterizations (See the 
last column of Table [2]). Since the RMF model we use in this work gives the same isoscalar 
and isovector effective masses, we set them equal in the reference SHF model too. Finally, the 
isoscalar parameters of the two reference Skyrme forces are then re-adjusted to fit the binding 
energy and charge radius of 208 Pb by adjusting only the saturation density po and the binding 
energy Eq of SNM. These ensure that the predictions for the charge radii and binding energies 
of other doubly closed-shell nuclides will be within 1-2% accuracy. In terms of the predicted 
values of isoscalar and isovector bulk observables, both corresponding RMF and SHF models 
are therefore almost equivalent. 

Having obtained the PNM optimized parameter sets, the la errors on these two purely 
isovector parameters can be translated into equivalent errors on the symmetry energy parameters 
and the neutron skin thickness of 208 Pb using the covariance analysis (See Table [3]). The errors 
in J are less than ±1 MeV for all the parameterizations. The RMF model gives a relatively 
small error in L of around ±2 MeV, while the SHF model gives a much larger error around ±6 
MeV. Table [3] appears to indicate that within the la errors, both models are consistent in their 
predicted values of J and L. However, a la joint confidence regions in the J-L plane plotted 
in Fig. EJa) for both RMF and SHF models shows that in fact the two models predict non- 
overlapping regions in J-L space. Both models show a positive correlation between J and L, but 
with differing slopes. The origin of this difference lies mainly in the values of the higher-order 
symmetry energy parameters that are predicted upon optimization. There is a strong model 
dependency in the prediction for the curvature parameter of the symmetry energy K sym (see 
Table ED. For example, after the PNM optimization IU-FSU predicts K sym = -6.8 ± 12.9 MeV, 
while its Skyrme-like version predicts a smaller value of K sjm = —130.2 ± 13.3 MeV. When we 
plot the la joint confidence regions in the K sym -L plane for both RMF and SHF models [see 
Fig- El (b)] further differences can be seen: there is, generically, a negative correlation between 
the slope of the symmetry energy and K sym in the RMF model, while this correlation is positive 
in the case of the SHF model. 



Table 3. Isovector observables and associated la error bars from four reference 
parameterizations after the PNM constraints are applied. Values are shown for the symmetry 
energy at p = 0.1 fm -3 50.1 and at saturation density J, slope parameter L, curvature parameter 
-Ksymj isospin-dependent part of incompressibility K T , and the neutron skin thickness i? s kin of 



208 



Pb. All the quantities are in MeV apart from i? s kini which is in fm. 





5*0.1 


J 


L 


K 


K T 


-Rskin 


NL3* 
SkNL3* 
IU-FSU 
SklU-FSU 


24.9 ± 0.4 
24.5 ± 0.3 
24.9 ± 0.4 
24.4 ± 0.3 


30.7 ± 0.7 
31.0 ± 0.9 
31.4 ± 0.7 
31.4 ± 0.9 


50.3 ± 1.8 

46.4 ± 6.4 
52.9 ± 2.0 
48.0 ± 6.2 


39.2 ± 17.8 
62.7 ± 16.6 
-6.8 ± 12.9 
-130.2 ± 13.3 


-284.6 ± 29.4 
-380.0 ± 15.2 
-257.6 ± 22.3 
-343.9 ± 15.3 


0.18 ± 0.01 
0.16 ± 0.01 
0.18 ± 0.01 
0.16 ± 0.01 



56 - 



52 - 



CD 



48 - 



44 - 



40 - 



IU-FSU 




(a) 



IU-FSU^ 



_L 



'SklU-FSU 



_L 



(b) 



_L 



_L 



- 56 



- 52 



- 48 



- 44 



- 40 



CD 



30 31 32 33-150 -100 -50 
J (MeV) K um (MeV) 



sym 



Figure 2. la joint confidence regions for the symmetry energy J and its slope parameter L 
(a), and for the slope parameter L and curvature parameter K sym (b) of the symmetry energy 
at saturation density from the IU-FSU and SklU-FSU parameterizations. 



Particularly attractive are predictions for the neutron skin thickness of Pb after the 
PNM optimization. The original NL3* parametrization predicts a very thick neutron skin of 
i?skin = 0.29 fm, while its Skyrme counterpart predicts a value of -R s kin = 0-27 fm. On the other 
hand, both the original IU-FSU and SklU-FSU parameterizations predict a much lower value of 
-Rskin = 0.16 fm for 208 Pb. These values are consistent with the current experimental result of 
-Rskin = 0.33lg'Jg fm for the neutron skin thickness of lead obtained using electroweak probes in 
the PREX experiment [37]. After the PNM optimization, we find that in general the RMF model 



predicts i? s kin = 0.18 ± 0.01 fm for i08 Pb, while the SHF model predicts i? skin = 0.16 ± 0.01 fm 
(See Table [3]), which are marginally consistent with each other within the la error-bars. Thus 
both models generically predict a thin neutron skin thickness of lead. Moreover, the error-bars 
coming from this constraint are very small. This is a particularly provocative result, since if the 
new PREX experiment reduces the error bars without moving the central value for the neutron 
skin, almost all current models of the nuclear structure would need to be modified. Also, this 
would appear to call for a significant modification of the PNM microscopic calculations. 




Finally, we apply our results to the neutron star structure, by extrapolating the EoS from 
the RMF and the SHF models to higher densities. The neutron star matter is assumed to be 
charge neutral and in the /3-equilibrium condition with neutrons, protons, electrons, and muons. 
No exotic degrees of freedom are assumed. The equations of state from the four parametrization 
are then utilized in the general relativistic equation of stellar structure (known as the Tolman- 
Oppenheimer-Volkoff equation) to obtain the mass-vs-radius relation. In Fig. [3] (a) we display 
mass-vs-radius relations as predicted by the four original RMF and SHF parametrizations that 
predict a wide range of results for low mass neutron star radii. This can be mainly attributed 
to the density dependence of the symmetry energy, which is quite different in the original two 
parameterizations. Once calibrated to the PNM results, this difference almost vanishes within 
the same model as shown in the Fig. [3](b), i.e., both RMF and SHF parameterizations now match 
each other more closely. The only difference at high masses between the RMF parametrizations 
is now due to the £ parameter that results in a stiffer EoS of SNM in NL3* parameterizations 
at several times saturation density. Although both NL3* and IU-FSU parameterizations in a 
given RMF or SHF model predict similar radii, there is a clear difference between the RMF and 




Figure 4. Density dependence of symmetry energy from the four parameterizations after the 
PNM optimization (Figure is taken from Ref. PQ). 



the SHF predictions as a whole. In the case of IU-FSU and SklU-FSU we have almost a ~ 1 
km difference for the radius of a canonical neutron star. This discrepancy is even larger in the 
case of NL3*, which is about ~ 1.8 km. Thus, there is a strong model dependence when the two 
models are applied to neutron star structure calculations after the same PNM optimization. This 
model dependence primarily comes from the different density dependence of symmetry energy 
at supra-saturation densities. Looking at Fig. [J] we see that the symmetry energy is almost the 
same in all the models up to ~ l-5po. The RMF model predicts a monotonic increasing function 
of density for the symmetry energy, while the SHF model tends to give a decreasing symmetry 
energy at higher densities. We thus show that the low-density PNM constraints alone result in 
a distinct model dependency of radius predictions. With the similar saturation properties of 
nuclear matter constrained by the PNM EoS, one can obtain different radii for a given neutron 
star mass. Although the PNM optimization tightly constrains the symmetry energy up to a 
little above the saturation density, its high density behavior that is very crucial in determining 
neutron star radii, still remain unclear. To elucidate this long-standing problem further we need 
to rely on the heavy-ion collision experiments [38[ 1391 HQ] and neutron star observations [411 [4"2] . 

6. Summary 

Using our best knowledge of the PNM EoS below and around saturation density from ab initio 
calculations, we constrain the density dependence of the symmetry energy for the RMF and 
SHF models, by optimizing the two pure isovector parameters from each model while keeping 
the values of other parameters so that the errors of predicted binding energies and charge radii 
of medium to heavy nuclei remain to be less than 2%. 

We show that such fits result in very similar predictions for the symmetry energy J and its 
slope parameter L at saturation density from both models as long as the nucleon effective mass 
from both RMF and SHF models is chosen to be consistent [36, 35J. When the error bounds 



are plotted as ellipses in the J-L plane, a positively-correlated relationship between J and L 
is observed for both models. However, different slopes are obtained from the RMF and SHF 
models, and the two ellipses have no overlapping area in the plane. This model dependence 
comes from the different values of K sym and higher-order symmetry energy parameters. 

Predictions of a neutron skin thickness -R s kin for 208 Pb are similar from both models and are 
within the error-bar of the latest experimental data. Both models predict a very thin neutron 
skin. Although the PNM constraints lead to broadly similar behaviors of the symmetry energy 
as a function of density up to « 1.5po> they deviate significantly at higher densities due to the 
differences in the functional form of the symmetry energy. This results in a striking difference 
in the predictions of the neutron star radii from both models. 

We showed the possible differences in predictions from two phenomenological energy-density 
functional forms, when the same experimental or theoretical constraints up to saturation density 
are applied. Care must be taken interpreting observational and experimental constraints from 
different nuclear models, and searching for a robust and better-determined EDF is necessary. 
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